clear; clc; close all;

u   = 0.1;   % 非线性控制参数
k   = 1.8;   % 耦合强度
a   = -1;    % 忆阻器参数 a
x0  = 0.1;   % 初始值 x0
q0  = 0.1;   % 初始值 q0

Ntotal   = 1200;  % 总迭代步数
Ndiscard = 200;   % 舍弃前 200 步
Nkeep    = 100;   % 分岔图只取最后 100 步

% 扫描 b 从 0 到 40
bVals = linspace(0, 40, 601);

% 用于存放李指数
LE1 = zeros(size(bVals));
LE2 = zeros(size(bVals));

% 用于存放分岔图数据 (b, x)
bifData = [];

for i = 1:length(bVals)
    b = bVals(i);

    % 1) 调用 LEs2D 计算李雅普诺夫指数
    [le1, le2] = LEs2D(@nhmm, u, k, a, b, x0, q0, Ntotal, Ndiscard);
    LE1(i) = le1;
    LE2(i) = le2;

    % 2) 迭代系统, 只取后 Nkeep 步的 x
    [xArr, ~] = iterate_nhmm(u, k, a, b, x0, q0, Ntotal);
    xSegment  = xArr(end - Nkeep + 1 : end);

    % 收集 (b, x) 用于分岔图
    bifData = [bifData; b * ones(Nkeep,1), xSegment(:)];
end

% (a) 李指数 vs. b
subplot(2,1,1);
plot(bVals, LE1, 'r-', 'LineWidth',1); hold on;
plot(bVals, LE2, 'b-','LineWidth',1);
yline(0,'--k','LineWidth',1);
xlabel('b');
ylabel(' LEs ');
title(' LEs 指数随 b 的变化');
grid on; box on;
legend('\lambda_1','\lambda_2','Location','best');

% (b) 分岔图 (x vs. b)
subplot(2,1,2);
plot(bifData(:,1), bifData(:,2), '.', 'MarkerSize',2, 'Color',[0.5,0,0.5]);
xlabel('b');
ylabel('x');
title('分岔图');
grid on; box on;
